@roystgnr Please take a look to see if the update makes sense. It does solve the issue I am having now with MOOSE but I am not sure if it would affect other stuffs. Thanks.

Job Coverage on bfa2615 wanted to post the following:

Coverage

Diff coverage report

Full coverage report

This comment will be updated on new commits.

Job Test MOOSE debug on 30c1237 : invalidated by @miaoyinb

Definitely makes sense. Using an absolute tolerance rather than a relative one bugs me (at least in that first test; t is non-dimensionalized already so "absolute" is already relative), but we can scale by the ray or edge length to fix that.

What worries me more is: IIRC, that (t >= 0 && t < 1) test was very carefully chosen, with the t==0 case needing special handling because of the "intersects at exactly one point" case, and turning that from an exact test into a fuzzy test might have broken some of those cases. It's at least a very good sign that we're still passing all unit tests ... but the very fact that you're running into problems means that there are still corner cases which our unit tests don't cover.

Can you give an example (in the form of a new unit test would be ideal, but if you just have a MOOSE input file with XYDelaunay then I'll work from there) of a case that fails with the current code but works with this? That would be helpful in figuring out if we can improve the solution here (in particular, I could try to rescale that absolute tolerance, but I have no idea whether the result would still fix your failing case!) or if I can find a better solution. And a new test case would be indispensable in making sure we never revert the fix. Imagine if someone later does run into another corner case, and their changes effectively revert this in order to fix it; I'd have no way to realize (except maybe a twinge of deja vu when I say to them "It's at least a very good sign that we're still passing all unit tests") that we were breaking you again.

Sure. I am developing a complex mesh generator that uses XYDelaunay as sub-generator and met with this issue. Let me see if I can make a simpler example.

This is a simple input that has the issue with the current code but can be run without issue after the modifications here.
I did notice that the issue is caused by a hole boundary node at {x1,0,0} and an external boundary node at {x2,0,0}, because a {1,0,0} ray direction is used here.

[Mesh]
  [pccmg1]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 5
    num_sectors_per_side = '4 4 4 4 4'
    background_intervals = 2
    polygon_size = 8.0

    ring_radii = '6.0 7.0'
    ring_intervals = '2 2'
  []
  [bdg1]
    type = BlockDeletionGenerator
    input = pccmg1
    block = 4
    new_boundary = 10000
  []
  [tg1]
    type = TransformGenerator
    input = bdg1
    transform = ROTATE
    vector_value = '90.0 0.0 0.0'
  []  
  [pcg]
    type = ParsedCurveGenerator
    x_formula = 'R*cos(t*2*pi)'
    y_formula = 'R*sin(t*2*pi)'
    section_bounding_t_values = '0 1'
    nums_segments = 60
    constant_names = 'pi R'
    constant_expressions = '${fparse pi} 100'
    is_closed_loop = true
  []
  [xyd]
    type = XYDelaunayGenerator
    boundary = pcg
    holes = tg1
    refine_boundary = false
    stitch_holes = 'true'
    desired_area = 10
  []  
[]

That's perfect; thanks. I can replicate easily with that.

I'm still hoping to come up with a better fix than your patch here, but I'm also still failing to. If I don't come up with something soon, we can just merge your patch but add a unit test based on your failing case here to make sure future revisions don't break it again. Mind if I push to your branch?

Sounds good. Please feel free to push to my branch.

@roystgnr I think your concerns make sense. I did meet with a false error even with this modification. I have not yet been able to replicate the new error using a simple input. But I think it should be related to the situation you mentioned.

I am thinking if we should use a random ray direction instead of {1,0,0}; or use two ray directions; or do a prescreening of ray direction to make sure the ray is away from a node. I was just thinking out loud here.

Well, I had been hoping this was reliable enough to merge, even if it wasn't a platonic ideal yet, so long as we had a unit test covering the previously-broken case ready for the next time someone wants to take a shot at that ideal. But I hadn't seen your latest comment yet.

I'd be afraid of a random or a "best out of 3" or whatever ray search just covering up the problem, making it harder to replicate, rather than actually fixing it. Nothing but p^n odds stops any other rays from hitting the same sort of bug.

For your new failure case, if you can come up with an input file that triggers it with standard MOOSE, even if it's not simple, share it and I'll take a crack at simplifying it as well?

The best I can do so far is to provide an input using the existing MOOSE MGs and load the mesh created by my under-development MGs from exodus files here (accs.zip).

The example runs well with the original libMesh but fails the hole verification with the modified libMesh.

assm_pitch = 100
assm_side = '${fparse assm_pitch/sqrt(3.0)}'
[Mesh]
  [acc1]
    type = FileMeshGenerator
    file = acc1.e
  []
  [acc2]
    type = FileMeshGenerator
    file = acc2.e
  []
  [hex]
    type = PolyLineMeshGenerator
    points = '${fparse assm_side/2.0} ${fparse assm_pitch/2.0} 0.0
                  ${fparse assm_side} 0.0 0.0
                  ${fparse assm_side/2.0} ${fparse -assm_pitch/2.0} 0.0
                  ${fparse -assm_side/2.0} ${fparse -assm_pitch/2.0} 0.0
                  ${fparse -assm_side} 0.0 0.0
                  ${fparse -assm_side/2.0} ${fparse assm_pitch/2.0} 0.0'
    loop = true
    num_edges_between_points = 10
  []
  [xyd]
    type = XYDelaunayGenerator
    boundary = hex
    holes = 'acc1 acc2'
    refine_boundary = false
    stitch_holes = 'true true'
    desired_area = 10
    # verify_holes = 'false'
  []
[]

Well, I can replicate this one too, so that's a good start!

And it's exactly what I was fearing might be a failure case: our ray from one hole point is "intersecting" (within tolerance) at one vertex of the other hole, but this time it's a "tangent" intersection rather than a "through" intersection, and our "check for double-intersections and don't count them" code isn't working after your patch.

I'll turn it into a unit test and we'll see what we can do.

Got some time to re-investigate this. My idea is to add a random, small but non-zero shift to the ray direction. It seems to work for my cases but not sure if it is generally a good idea. (obviously it failed some tests 😢) Any insight would be greatly appreciated.

I think a random shift is a terrifying idea - it won't actually fix the bug, it'll just make the bug harder to replicate.

My vague intuition here is that we've hit a problem where using any tolerances is a bad idea. The goal isn't to catch corner cases here and there, it's just to treat corner cases consistently. Whether we think we're hitting a boundary vertex or not doesn't actually matter. What matters is that, when we look at the same boundary vertex from each of the edges sharing it, we always come to the exact same conclusion in both cases about whether our ray intersects it or goes "above" it or "below" it.

Hmm... maybe that's what I'm doing wrong here? Instead of trying to find the parameter u for the distance along each of the two edges our intersection is at, then check for u∈[0,1], maybe we should be testing whether a ray is clockwise/counterclockwise/intersecting each of the three vertices of those two edges. By doing a vertex test instead of an edge test we'd guarantee consistency.

I think a random shift is a terrifying idea - it won't actually fix the bug, it'll just make the bug harder to replicate.

My vague intuition here is that we've hit a problem where using any tolerances is a bad idea. The goal isn't to catch corner cases here and there, it's just to treat corner cases consistently. Whether we think we're hitting a boundary vertex or not doesn't actually matter. What matters is that, when we look at the same boundary vertex from each of the edges sharing it, we always come to the exact same conclusion in both cases about whether our ray intersects it or goes "above" it or "below" it.

Hmm... maybe that's what I'm doing wrong here? Instead of trying to find the parameter u for the distance along each of the two edges our intersection is at, then check for u∈[0,1], maybe we should be testing whether a ray is clockwise/counterclockwise/intersecting each of the three vertices of those two edges. By doing a vertex test instead of an edge test we'd guarantee consistency.

Makes sense. The random shift would make the odds of hitting this bug nearly zero but not actually zero...

I agree we should try a different approach other than edge test.